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Abstract . 


This paper describes a finite element formulation 
designed to simulate polymer melt flows in which both con- 
ductive and convective heat transfer may be important, and 
illustrates the numerical model by means of computer experi- 
ments using extruder drag flow and entry flow as trial prob- 
lems. Fluid incompressibility is enforced by a penalty 
treatment of the element pressures, and the thermal convec- 
tive transport is modeled by conventional Galerkin and opti- 
mal upwind treatments. 

Introduction 


Process development for polymer melt flow operations has 
in the past been largely empirical, due in part to the dif- 
ficulty of obtaining closed-form solutions to the governing 
equations which satisfy the complicated boundary conditions 
encountered in actual equipment. The finite element method 
promises a significant advance in our capability to obtain 
numerical solutions to these processing problems, in that it 
is able to generate approximate solutions for irregularly 
shaped boundaries. In an earlier paper (_1_) this author 
described the utility of "penalty" finite elements in pre- 
dicting the velocity and stress fields developed in melt 
processing, and this paper will describe the extension of 
this method to include coupled heat transport effects. It 
is hoped that such a capability will help the process 
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designer avoid such processing problems as thermomechanical 
degradation or improperly sized equipment, and further to 
process materials so as to produce optimal microstructures 
and properties. 

Formulation 


The analyses to be reported here were obtained by means 
of isoparametric penalty elements developed at MIT and 
incorporated within the finite element analysis program FEAP 
listed in the text by Zienkiewicz (2^. The element formula- 
tion will be described briefly here, and the reader new to 
finite element methodology is referred to standard texts orr 
the subject (2_,_3) for a more extensive discussion of the 
various concepts. The description here will be that for 
plane flow problems, but the concepts are readily extended 
to axisymmetric and three-dimensional formulations- as well. 

The solution domain of the problem is discretized as an 
assemblage of elements along whose boundaries are located 
nodes at which velocities and temperatures T^ are sought. 
(The i subscript ranges over the element nodes, and under- 
lining indicates a vector or matrix quantity.) The nodal 
values of the x and y components of velocity and the temper- 
ature are listed in a vector _a^ : 

Si = K.Vi.Ti) = (Mi.Ti) 
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The values of velocity and temperature at an arbitrary 
position within an element can be approximated by interpo- 
lating among the nodal values: 

u(x,y) = N i ( x ,y) u ± 

. T( x , y ) = N i ( x ,y) T ± 

The interpolation functions N^(x,y) are obtained from stand- 
ardized subroutines and can be selected to provide linear, 
quadratic, or higher-order interpolation among the nodal 
values. This interpolation is a central concept in finite 
element analyses, in that the actual variation of the 
unknowns within an element is replaced by a low-order poly- 
nomial interpolation. 

The rate-of-deformation tensor V can be obtained from 
the velocity field as: 

V T = (3u/3x, 3v/3y, 3v/3x + 3u/3y) 

V = Lu = LN i u i = B i u i 



3/3x 

0 


N i,x 0 

L = 

0 

3/3y 

®i = 

0 N i ,y 


3/3y 

3/3 x 


- N i ,y N i,x. 


In the case of Newtonian fluids of viscosity y, the devia- 
toric stresses t can be obtained from V as: 


4 



( T T T ) 

xx » yy » xy' 
t = D V 


D = y 


2 0 0 
0 2 0 
L0 0 1 


The hydrostatic stress p can be obtained by considering the 
volumetric dilation V y : 

V v = V T u = m T B i u i ' 

V T = (3/3x, 3/ 3y ) , m T = (1 , 1,0) 

P = XV v = Xm T B i u i 

Here X is comparable to the fluid bulk modulus, and in the 
case of incompressible fluids can be regarded as a parameter 
which "penalizes" fluid compressibility. By appropriate 
choice of X (generally X is chosen as 10^y (4_)), incompres- 
sibility can be enforced to whatever degree is desired. The 
total stress a, is then: 


a = t _ + mp 


The nodal unknowns _a^ are to be chosen so as to satisfy 
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the governing equations in an integral sense; this can be 
done by using a Galerkin weighted residual formulation of 
the conservation equations for momentum and energy trans- 
port: 



f 


(V T kVT 


Q) dft 


0 


Here ft is the element area, Q is the internal heat genera- 
tion, k is the thermal conductivity, and the weighting func- 
tions are the same functions used for interpolation. 
Using Green's theorem on the above integrals and substitut- 
ing the interpolated expressions for a and T, one obtains 
the matrix equation: 


Ka = f 


where 


K X + K y 0 


K = 


0 






(m T Bj) T A(m T B) dft 


K 


jk 




dft 
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K jk = /<’ TN j>KVN k 
f u = / Hl t* dr 

fT = J N i Q dft 


dfi 


The jk subscripts identify the influence of the velocity (or 
temperature) of the kth node on the force (or thermal flux) 
at the jth node, and the _t refers to tractions applied to 
the element boundary r. In these analyses Q is taken as 
arising only from viscous dissipation: 

Q = t T V 

The FEAP code is configured to cycle through each ele- 
ment in turn, evaluating the various K's and f's by Gauss- 
Legendre numerical integration and adding each element's 
stiffness or loading to the global arrays relating all of 
the problem's nodal unknowns (the a's) to the loading terms 
(the f's) . In the analyses to be presented below, four-node 
linear or eight-node parabolic quadrilateral elements were 
used in conjunction with 2x2 numerical integration on the K U 
and K terms. "Selective reduced integration" was used to 
evaluate the terms; this reduced integration order is 
necessary to prevent the large X parameter from driving the 
a vector to zero (2,4) . 

Iteration must be used if the coefficients in the formu- 
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lation (the y, k, Q) depend on the final solution for a. In 
the solutions to be given below, the y and k are taken as 
constant, and the velocity and temperature fields are cou- 
pled only through the Q term. Here only two iterations are 
necessary: the velocities and viscous dissipation are com- 

puted correctly in the first iteration, and the thermal val- 
ues in the second . 

Although momentum convection is generally negligible for 
polymer melt flows (Reynolds' numbers on the order of 10 - ^) } 
thermal convection may be significant. The ratio of convec- 
tive to conductive transport is given approximately by the 
Peclet number Pe : 

Pe = Uhpc p /k. 

Here U and h are a characteristic velocity and length, while 
p and c p are the fluid density and specific heat. Using 
typical values for down-channel extruder drag flow of low 
density polyethylene (5^), for instance, one computes Pe - 
5000 . 

Thermal convection can be incorporated in the finite 
element analysis by adding to the K T terms the quantity 

K jk = f N jP c pH T V N k dQ 

Incorporation of this term is straightforward, since the 
velocities u are known after the first iteration. The con- 
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vective terms are unsymmetric , however, which requires using 
less efficient storage and solution schemes than are possi- 
ble for symmetric systems. Even more seriously, inclusion 
of convective transport effects tends to produce numerical 
instability in the computed results, as will be illustrated 
below . 

Down-channel extruder flow . 

The capabilities of the coupled flow-heat transfer model 
will be illustrated first by means of a simulation of down- 
channel extruder flow, which may be modeled by the usual 
unwound-channel approximation as a straight rectangular 
channel containing fluid which is being dragged downstream 
by a plate covering the channel top (_1_,_5) . This is the 
well-known Couette flow problem, and is useful as a model 
problem both by virtue of its relevance to extruder flow and 
by the availability of theoretical solutions against which 
the numerical results may be checked. Figure 1 shows the 
normalized velocity and temperature distributions across the 
channel height, where the drawn curves are the theoretical 
predictions and the open symbols are the numerical values 
obtained using a 4x4 mesh of eight-noded parabolic elements. 
The theoretical temperature variation is given as 

T* = (Br/2)y* (1-y*) 
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Here y is the height from the root of the channel 

£ 

normalized on the total channel height, T = (T-T w )/t is 
the temperature relative to the wall temperature, and Br = 
UV 2 /k is a Brinkman number which provides an indicator of 
the importance of viscous dissipation relative to heat con- 
duction. The numerical results of figure 1 were obtained 
with Br = 8 . 

Figure 2 displays the results of a simulation of a 
slightly more complicated extension of the problem of figure 
1. Here the root and barrel temperatures are different, and 
the flow is retarded by the presence of a pressure gradient 
which is positive in the downstream direction. This simula- 
tion is more typical of extruder operation, in which the 
barrel but not the screw is heated and the constriction at 
the die produces a pressure gradient. 

Entry flow . 

Entry flow is a much-studied problem which provides a 
more demanding test of numerical codes than does the viscom- 
etric Couette flow problem described above. Theoretical and 
numerical aspects of this flow have been reviewed recently 
by Boger (6^), and comparison with these other results pro- 
vide a check of at least the flow aspects of the model. A 
reservior-to-capillary size ratio of 4:1 was selected, since 
the flow results become independent of size ratio at this 
value or larger. Figure 3 shows the grid of 100 four-node 
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linear elements used to model the upper symmetric half of 
this problem (273 degrees of freedom) , and the portion of 
the velocity field near the throat. A fully-developed 
Poiseuille velocity distribution was imposed on the entrance 
to the reservoir. The velocity distribution in the capil- 
lary becomes stable within a few diameters of the entrance, 
but the field is perturbed far upstream. Figure 4 shows the 
streamlines developed from the flow field by contour inte- 
gration along vertical lines through the grid, and normal- 
ized so as to be zero along the boundaries and unity along 
the centerline. These streamlines are identical with pub- 
lished experimental and numerical results, although the grid 
used here was not intended to be fine enough to capture the 
weak recirculation which develops in the stagnant corner of 
the reservoir. 

Figure 5 shows the contours of shear stress near the 
throat, normalized on the theoretical value expected at the 
capillary wall. The shear contours are in agreement with 
the isochromatic fringes seen in birefringence photography 
of entry flow (7), which is an example of a possible experi- 
mental verification of the flow model. 

Although certain grids using the penalty formulation are 
known to develop a pathology known as "checkerboarding" in 
which the computed pressures oscillate from element to ele- 
ment (40, this grid is free of that troublesome (but fixa- 
ble) result. Figure 6 .shows the profile of computed pres- 
sures along the elements adjacent to the centerline, 



normalized on the value expected in the reservoir in the 
absence of entry or exit pressure . losses . As is well known, 
the development of the flow field at the entrance to the 
capillary produces an additional pressure drop which must be 
accounted for in capillary viscometry by means of a "capil- 
lary correction factor". The correction factor, defined as 
he entrance pressure drop AP ent divided by twice the capil- 
lary wall shear stress, is given by Boger as 0.589 based on 
the numerical and experimental work of several workers. The 
value computed from the data of figure 6 is 0.594, which is 
within 1 % of Boger's value. 

As described above, the temperature field is computed 
using the finite element formulation of the heat conduction 
equation, with the viscous heat generation being computed 
from the stress and velocity fields obtained during the 
first iteration of the problem. The temperature contours, 

p 

normalized on the maximum centerline temperature T c = UV /3k 
expected for capillary Pois'euille flow, are shown in figure 
7. A hot spot is noted just upstream of the entrance, due 
to the combination of increased dissipation and greater dis- 
tance from the cooler boundaries at this position. This hot 
spot leads to higher temperatures both upstream and down- 
stream than would have been expected in simple Poiseuille 
flow, due to conductive heat transfer from the hot region. 

The temperature map of figure 7 is for Pe = 0, the case 
for which convective terms are not included in the heat 
transport calculations. Whereas all convective terms vanish 
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identically in the Couette flow described in the previous 
section, convective terms are present in the entry flow 
problem, and these might be expected to be significant for 
many real problems in the flow of polymer melts through 
dies. Figure 8 shows the profile of temperatures along the 
centerline for several values of Peclet number (Pe based on 
the capillary half-height and the maximum velocity) . It is 
seen that increasing convective transport serves to carry 
the cooler upstream particles downstream so as to diminish 
both the magnitude of the hot spot and its ability to con- 
duct heat upstream. As Pe increases, the influence of the 
entrance region is compressed into a smaller and smaller 
boundary layer at the throat, leading to numerical difficul- 
ties associated with the inability of the fixed mesh to cap- 
ture effects which occur in the smaller region. At Pe = 10, 
the values computed by the conventional Galerkin process 
(solid lines) suffer from noticeable instability, and by 
Pe = 100 the Galerkin solution has degenerated completely to 
meaningless values. 

The instability associated with convection analyses has 
plagued both finite difference and finite element analyses, 
and in recent years several ad hoc treatments have been pro- 
posed for obtaining stable and nonosc illatory results in the 
presence of strong convective terms. One of the most popu- 
lar of these treatments is "upwinding", in which the 
upstream portion of the interpolation is weighted more heav- 
ily in an attempt to smooth out the high local downstream 
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gradients which cause the problem in convective flows. The 
dotted lines in figure 8 were obtained using a convenient 
’’optimal" upwinding formulation suggested by Hughes (4), in 
which the upwind weighting is obtained simply by relocating 
the integration point at which the element shape functions 
are evaluated. The position of the upwind point is deter- 
mined for each element based on its local Peclet number, so 
that optimum results may be obtained. Figure 8 shows that 
the upwind solution for Pe = 10 is considerably smoother 
than the Galerkin result, and that the results for Pe = 100 
in which the upstream temperatures are carried well into the 
capillary are captured in a smooth and reasonable manner. 
The reader is cautioned, however, that the use of these 
upwinding procedures is still controversial ,‘ and is directed 
to a provocative paper by Gresho and Lee (8) in which the 
drawbacks of this approach are detailed. 

Conclusion . 


This paper has described an analytical tool which can 
predict velocities, stresses, and temperatures in noniso- 
thermal flow situations of the sort encountered in many 
polymer melt processing operations. Such a model cannot be 
expected to replace the experience and intuition which pro- 
vide the basis for most process design today, but it is 
hoped that this inexpensive and easily implemented model can 
provide a means by which the process' designer's intuition 
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might be expanded. Properly used, it can be a valuable 
additional tool at the process designer's disposal.. 

Development of this model is continuing in our labora- 
tory, and among the aspects still under development are 
capabilities for transient flows, reactive fluids, free sur- 
faces, and wall slip. Although incorporation of fluid elas- 
ticity is desired due to its importance in many polymer melt 
flows, such a development has proven elusive to a number of 
well qualified groups in the past several years. At pres- 
ent, it seems prudent to let the theoretical aspects of 
elastic effects be developed further before attempting their 
incorporation in a general process model. 
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Figure 1 - Normalized velocity and temperature distrib- 
utions in Couette flow. 
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Figure 2 - Couette flow with imposed pressure gradient, 
and with different root and barrel temperatures. 
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Figure 4 


- Entry flow streamlines. 
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